//
// Created by denglibin on 2019/7/10 0010.
//
/********解多元一次方程组******/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
/**
 * 消元计算最后一个未知数的值
 * 示例:
 * ax + by = f
 * cx + dy = e
 * 消去x后:
 *  bcy - ddy = fc- ed -> (bc-da)y=(fc-ea) -> y = (fc-ea)/(bc-da)
 * @param arr 系数数组（二维数组作为参数必须指定列长度，所以用一维数组代替二位数组）arr2[i][j] = arr[i * (n+1) +j]
 * @param n  几元 n+1表示列长度
 */
static double  elimination(double * arr, int n){
    if(n == 1){
        if(arr[0] == (double)0){
            return 0;
        }
        return arr[1] / arr[0];
    }
    //消元后（n-1元）的系数（二维数组 ）
    double * temp = calloc( (n-1) * n , sizeof(double));
    double a = arr[0]; //0行0列 （用第一行和后面的所有行作减法 最终得到n-1行）
    //消元 原来的n行n+1列变为n-1行n列 所以从1开始
    for(int i = 1; i < n; i++){
       for(int j = 1; j <= n;j++){
           double b = arr[0 * (n+1) + j ];//第0行第j列
           double c = arr[i * (n+1) + 0 ]; //第i行第0列
           double d = arr[i * (n+1) + j];  //第i行第j列
           temp[(i-1) * (n) + (j-1)] =  b * c -  d * a; //消元后 第i-1行第j-1列的系数
       }
    }
    printf("--------------------elimination-------------------\n"); //打印消元后
    for (int k = 0; k < (n-1) * n; ++k) {
        if(k % n == 0){
            printf("\n");
        }
        printf("%f ", *(temp+k));

    }
    printf("\n----------------------------------------------------\n");
    if(n > 2){ //多余2个未知数 再次消元
        double  value  =  elimination(temp, n - 1);
        free(temp);
        return value;
    }
    double  value = temp[1] /temp[0];
    free(temp);
    return value;
}
/**
 * 解多元一次方程组
 * ax + by = c
 * cx + dy = e
 * 二元则数组的2行3列
 * 三元数组则是3行4列
 * ....
 *
 * @param arr 系数和常数项的二维数组
 * @param n  n 元
 */
    static void equation_set(double *  arr, int n, double * result ){

    //消元解得最后一个未知数的值
    double  value = elimination(arr, n);
    result[n-1] = value;
    //多余一个未知数
    if(n > 1){
        //新的二位数组（因为已经解出一个未知数了，剩余未知数n-1个）
        double * newArr = calloc( (n - 1) * n, sizeof(double));
        for(int i = 0; i < n - 1; i++){ //n-1元
            for(int j = 0; j < n; j++){ //n 列
                if(j < n - 1){
                    newArr[i *( n ) + j] = arr[i * (n+1) + j]; //未知数系数和原数组一致
                }else{
                    //本行最后一个数(常数项)  最后一个未知数（已解出为：value）与系数相乘移到等式右边
                    newArr[i *( n ) + j] = arr[i * (n+1) + j+1] - arr[i* (n+1) + j] * value;
                }
            }
        }
        equation_set(newArr, n-1, result); //递归
        free(newArr);
    }
}

int main8(){
    int N = 4;
    srand(time(0));
    double  * arr = calloc(N *(N+1), sizeof(double));
    for(int i = 0; i < N * (N+1); i++){
        if(i % (N+1) == 0){
            printf("\n");
        }
        arr[i] = rand() % 20;
        printf("%f\t", arr[i]);

    }
    printf("\n");
    double * result = calloc(N, sizeof(double));
    equation_set(arr, N, result);

    printf("result:\n"); //结果
    for(int i = 0; i < N ; i++){
        printf("%f\n", result[i]);
    }

    printf("validate result:\n");
    //验证结果是否正确：
    for(int i = 0; i < N ; i++){ //N行
        double  sum = 0; //该行等式左边的值
        for(int j = 0; j < N; j++){
            sum += arr[i * (N+1) + j] * result[j];
            printf("(%f) * (%f)", arr[i * (N+1) + j] , result[j]);
            if(j < N-1){
                printf(" + ");
            }
        }
        double  rightValue = arr[i * (N+1) + N]; //改行等式右边的值
        printf(" = %f----%f\n",sum,  rightValue); //相等则验证通过

    }

    free(arr);
    free(result);

    return 0;
}

